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In open system approaches with non-Markovian environments, the process of inserting an indi¬ 
vidual mode (denoted as “pseudomode”) into the bath or extracting it from the bath is widely 
employed. This procedure, however, is typically performed on basis of the spectral density (SD) and 
does not incorporate temperature. Here, we show how the - temperature-dependent - bath correla¬ 
tion function (BCF) transforms in such a process. We present analytic formulae for the transformed 
BCF and numerically study the differences between factorizing initial state and global thermal (cor¬ 
related) initial state of mode and bath, respectively. We find that in the regime of strong coupling of 
the mode to both system and bath, the differences in the BCFs give rise to pronounced differences 
in the dynamics of the system. 
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I. INTRODUCTION 

In many open quantum system approaches, the micro¬ 
scopic model underlying the environment consists of an 
infinite number of harmonic oscillators linearly coupled 
to system degrees of freedom m- The flexibility of this 
particular model is owed to the fact that both spectrum 
and frequency-dependent coupling of the environmental 
modes can be adjusted to reproduce features observed in 
experiments, e.g. to describe the effect of polar solvents 
on dyes US] or to treat vibrational modes of molecules 
Hill. The spectrum and the frequency-dependent cou¬ 
pling can be encoded in a quantity called spectral density 
(SD), which is a real function of frequency. 

The influence of the environment on the system is even¬ 
tually determined by the bath correlation function (BCF) 
[UEJISHSl, which quantifies temporal correlations of envi¬ 
ronmental degrees of freedom. Via the initial state of the 
environment, environment temperature is incorporated 
in the BCF. 

The partitioning of the total system into explicitly 
treated degrees of freedom and a rest which is treated 
by means of a BCF is not fixed in the first place, and 
different choices might be expedient from different (com¬ 
putational) points of view [5]. That is, in the simu¬ 
lation of molecules embedded in an environment, one 
might explicitly include one or more strongly coupled 
(important) vibrational modes in the system part and 
treat the remaining part effectively. Alternatively, one 
could consider only the electronic degrees of freedom of 
the molecules as system part, which in turn are coupled 
to a highly structured environment. The very same con¬ 
siderations apply in the case of some electronic degrees 
of freedom being coupled to an imperfect (lossy) cavity 
[SlIIDl. While an unstructured environment is more easily 
handled in simulations, the rapid growth of the Hilbert 
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space associated with the explicitly treated environmen¬ 
tal modes makes the second choice the favorable one. 

The effect of different system-environment partitioning 
has already been discussed in literature EUnHSI, how¬ 
ever, the discussion mostly focused on the SD. In this 
work, we study how the BCF transforms under differ¬ 
ent system-environment partitionings. In particular, we 
examine the effect of two different initial states, a factor¬ 
izing and a correlated one, on the transformed BCF. 

We consider an exemplary model consisting of a har¬ 
monic bath coupled to a single mode, called “pseudo¬ 
mode” (PM). (See Ref. [5|for a discussion of how a PM 
relates to a vibrational mode.) This PM we then cou¬ 
ple to another harmonic oscillator that acts as a system, 
thus allowing us to study the influence of the different 
BCFs on a system dynamics. For strong coupling be¬ 
tween PM and both system oscillator and ohmic bath we 
find pronounced differences in the dynamics of the mean 
occupation number of the system oscillator, thus stress¬ 
ing the importance to take heed of the initial state of the 
composite system in this regime. 

The paper is structured as follows: In Sec. |n| we in¬ 
troduce the microscopic model on which our discussion 
is based. We outline the procedure according to which 
BCFs transform and state analytic formulae for the case 
of a single PM coupled to a harmonic bath. In Sec. [ml 
we evaluate the transformed BCFs numerically and dis¬ 
cuss some examples, highlighting the regimes in which 
notable differences are induced by different initial states. 
Finally, we summarize our findings in Sec. EYI 

Details of calculations are given in three appendices: 
In Appendix we review the definition of BCF and 
SD on basis of the microscopic model introduced in the 
main text. In Appendix]^ we explain how the PM and 
bath operators can be transformed into a basis in which 
the combined Hamiltonian of PM and bath is diagonal. 
Lastly, we review in Appendix[^an alternative derivation 
of the transformed BCFs on grounds of the Heisenberg 
equations of motion and state a numerical recipe to solve 
the occurring integro-differential equation. 
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II. MODEL SYSTEM AND ANALYTIC 
TRANSFORMATIONS 


In this section, we detail model Hamiltonian and 
framework necessary to perform the analytic transfor¬ 
mation of the BCF presented at the end of the section. 
To that end, we first review the standard model of a sys¬ 
tem linearly coupled to an environment of independent 
oscillators (Sec. IIA) and introduce the model Hamilto¬ 


nian we consider in this work (Sec. HB|. Subsequently, 


we discuss two particular ways to partition this Hamilto¬ 
nian into a system and an environment part (Sec. |H C ) 
and explain how the transformed BCF can be calculated 
(Sec. [HP I. Finally, we specify two initial environment 
states and present the corresponding transformed BCFs 
(Sec.[HE|. 


A. General properties of a system linearly coupled 
to an environment of independent oscillators 

We start our discussion by reviewing the standard 
model of a system linearly coupled to an environment 
of independent oscillators, which will allow us later on 
to point out differences in the transformed BCFs most 
clearly. 

In the standard model, the total Hamiltonian is parti¬ 
tioned into three parts, 

^tot + "^s-E +'He, (1) 

where Us denotes the Hamiltonian of the system, con¬ 
taining the degrees of freedom in which one is interested 
in, and 


= ( 2 ) 

the Hamiltonian of the environment. Here, is the anni¬ 
hilation operator of an environment mode with frequency 
An environment of the form of "He we call diagonal, 
meaning that the oscillators are uncoupled. 

The Hamiltonian accounting for the interactions be¬ 
tween system and environment reads 


"Hs-e 




ct 

mV 


H.c. 


( 3 ) 


The relevant quantity typically entering open quan¬ 
tum system approaches like Redfield [DE], Caldeira- 
Leggett [IlITl] and Non-Markovian Quantum State Dif¬ 
fusion (NMQSD) [SI [TS] is the bath correlation function, 
— for Hermitian Ls — given by 

= TrE {(C(t) + C't(t)) (C(f') + C\t')) pe(0)} , 

( 5 ) 

with C{t) defined as 

C{t) = ^ {k;c,) ^ ^ k;c,{t). (6) 

To obtain Eq. © in the given form, the total initial state 
is taken to be 


Ptot(O) = /5s(0) O pe(0). (7) 


This implies that no correlations between system and en¬ 
vironment exist before the interaction between system 
and environment is ‘turned on’. This assumption, which 
is typically introduced by virtue of the ease of computa¬ 
tion, establishes the significance of system-environment 
partitioning. 

For a non-Hermitian system operator Ls, the BCF is 
no longer given by Eq. 0 . Rather, two correlation func¬ 
tions are required m, reading 

ai{t,t') = TrE {C{t)C\t')p-E{0)} and ( 8 a) 

a2{t,t') = TrE {Ct(t)C(F)pE( 0 )} . ( 8 b) 


If the BCF is stationary, i.e., if is a function 

of the time difference only, = a{t — t',0), it is 

convenient to write a{T) = a(r, 0), with t = t — t'. Note 
that the stationarity of the BCF depends on the initial 
state of the environment in general. 

If the initial state of the diagonal environment /5e(0), 
which enters Eq. ([^ , is a thermal state (and if Ls is Her¬ 
mitian), the BCE of the environment is of the ‘standard’ 
form 


if(r) = J dujJ{oj) (^coth cos{ujt) — 


(9) 


with T = t — t'. In Eq. T is the temperature of the 
environment and J{uj) the SD. Eor a detailed review of 
SD and BCF in the standard case, see Appendix [A| 


B. Model Hamiltonian with PM 


The coupling Hamiltonian Hs-e couples the environ¬ 
ment modes via Ls linearly to the system, with 
strength k^. It is convenient to introduce the so-called 
spectral density encoding this frequency-dependent cou¬ 
pling as (w > 0 ) 

(4) 

M 


We now introduce the total Hamiltonian which we fo¬ 
cus on in this work, which is of the form 


LLtot = "Hrel + 'Hrel-PM + 'HpM + 'HpM-B + LLb- (10) 


In Eq. (10), 'Hrei contains the relevant degrees of free¬ 
dom we are interested in. This relevant part is via the 
Hamiltonian 


Note that we set h= k^ = 1 throughout this work. 


Hrel—PM 


( 5 VLI -f H.c.) 


( 11 ) 
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linearly coupled to the PM, whose Hamiltonian reads 

^PM = (12) 

Here, L is some operator in the Hilbert space of the rele¬ 
vant Hamiltonian 'Hrei, 9 a coupling constant quantifying 
the strength of the coupling, and VL the frequency of the 
PM with annihilation operator b. In addition, the PM is 
coupled to a diagonal bath 

1-Lb = '^ujxa[ax, (13) 

A 

where ujx are the frequencies belonging to the bath modes 
A with annihilation operators ax. The coupling Hamilto¬ 
nian HpM-B is taken to be bilinear, 

'HpM-B = ^ {^xaxb^ + H.C.), (14) 

A 

with Kx being the coupling constants quantifying the cou¬ 
pling between PM and bath modes. The generalization 
of our discussion to several PMs is straightforward in 
many cases of interest (e.g., for single linear chains of 
PMs [miiH] and multiple linear chains m), and will be 
addressed at the end of the section. 


C. System-environment partitioning 

We now consider two particular examples of assigning 
the PM to the different parts of the total Hamiltonian, 
illustrated in Fig.j^ This leads to different choices of the 
system Hamiltonian "Hs, the environment Hamiltonian 
■He, and the coupling between them. We denote the two 
different ways of partitioning by SI and SH for the system 
and El and EH for the environment, respectively. 




FIG. 1. Illustration of two different ways of performing a 
system-environment partitioning in the presence of a PM lin¬ 
early coupled to both system and bath. In (I), the system part 
SI consists of the relevant system degrees of freedom linearly 
coupled to the PM, which in turn is coupled to an (unstruc¬ 
tured) environment El whereas in (II) the system SII directly 
couples to a (structured) environment Eli including the PM. 


while the environment EH contains both PM and bath, 

HeII = HpM + HpM-B + Hb- (19) 

Accordingly, the coupling between system and environ¬ 
ment is given by 

Hsii-eii = Hs-pm- (20 ) 

Note that the resulting Hamiltonian is not in the stan¬ 
dard form of Eq. Q , as the environment is not diagonal. 
It can, however, be diagonalized by a simple transforma¬ 
tion, as detailed in Appendix [B| 


1. PM in the system 


D. Calculation of the BCF 


The first partitioning [(I) in Eig. is to take the PM 
as part of the system, which amounts to setting 


and 


Hsi = Hrel + Hrel-PM + HpM, 

Hei = Hb, 


’^SI-EI — HpM-B- 


(15) 

(16) 

(17) 


Note that the system now contains beside the relevant 
degrees of freedom also the PM, and that the environ¬ 
ment is in the standard form of Eq. Q. 


2. PM in the environment 


The second partitioning is illustrated in panel (H) in 
Fig.a Here, the system is given by 

(18) 


In terms of the SD, it is known how to transform be¬ 
tween different ways of system-environment partitioning 
0 [3 EH m [in]- Eor this procedure, it is sufficient to 
know the total Hamiltonian, as the SD is fully encoded 
in Htot- 

The BCE, however, depends on the environment state, 
denoted by /3ei(0) and /5eii(0) respectively for the two 
settings (I) and (H) in Fig. ^ Note that to obtain a total 
initial state of form Eq. ([7), we need /3tot(0) = /5si(0) ® 
/3ei(0) in setting (I) whereas ptot(O) = psii(O) (8)pEii(0) 
in setting (H). 

Microscopically, the BCE is the (two-time) correlation 
function of the environment operators in the system- 
environment coupling. In setting (I), the environment 
is diagonal and we can therefore directly use Eq. (§. For 
setting (H), in contrast, the environment EH is not diag¬ 
onal. Nonetheless, we can similarly to Eq. write down 
the time evolution of the environment coupling operator. 


B{t) = g* e 


H-Le 


be 


— z'HeII* 


= g*b{t), 


Hsil = Hrel 


( 21 ) 
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whose time dependence arises via transformation into the 
interaction picture with respect to 'Heii- 

For a Hermitian system operator, L = L\ and 'Hs-pm 
can be written as 'Hs-pm = L{g*b + which has a 
Hermitian environment coupling operator. In this case, 
the BCF is given by 

a(t, t') = TrEii { {B{t) + B\t)) {Bit') + B\t')) /3eii(0)} 
= {{B{t) + B^{t)) {B{t') + B^{t')))^,, , (22) 


This particular initial state is widely used, owing to its 
convenient properties in analytic calculations. The phys¬ 
ical assumption implied is that during thermal equilibra¬ 
tion with an ambient heat bath no correlations are built 
up between PM and bath El; a reasoning that only holds 
in the limit of vanishing coupling between PM and bath 
since only in this limit independent thermal equilibration 
(i.e., equilibration to the respective canonical states) of 
two coupled systems exists (cf. also Refs. [2S] and 


where peii(O) denotes the initial density operator of the 
environment and the subscript Eli of the trace indicates 
that the trace is taken over the environmental degrees of 
freedom. 

To evaluate the BCE ( [2^ , it is convenient to take ad¬ 
vantage of the existence of a linear transformation be¬ 
tween the PM operator b and the operators in which the 
Hamiltonian 'Heii and the initial state, respectively, are 
diagonal (cf. Appendix]^ . Specifically, we first transform 
the PM operator into the basis in which the environment 
Hamiltonian Heii is diagonal, b{t) = [S'c(t)]o, by means 
of the transformation matrix S. The time evolution of 
the annihilation operators of the diagonal Hamilto¬ 
nian Heii, however, is simply given by c^{t) = 
Subsequently, the operators are transformed into the 
basis in which the initial state is diagonal, if necessary, 
and the BCE is evaluated. 


Evaluating the BCF Eq. 
state, we find (see Appendix 


(22) with factorizing initial 
^for details) 






.{ijJri) 




+ (n(a;,) + I) 


(24) 


Here, n(w) is the mean occupation number of an envi¬ 
ronment oscillator with frequency w. 


n(a;) 


I 

(e/3<^ - I) ■ 


(25) 


Note that Eq. (24) is not in the form of Eq. ([^; in fact, 
we cannot even write a(t, t') = a{t — t', 0). 


E. Choice of initial states of the environment 


2. Thermal (correlated) state of PM and El 


As discussed in the previous section, the total ini¬ 
tial state in case (I) is typically taken to be /5tot(0) = 
/3si(0) 0 Pei(O). When moving the PM from the system 
part to the environment part, i.e., going from (I) to (H), 
one could thus reason that the initial state of the environ¬ 
ment EH should be given by /5eii(0) = /5pm(0) 0 /5ei(0). 
Conversely, if one considers the PM to be part of the 
environment EH from the very beginning on, there is no 
reason why the PM should be uncorrelated with EL Erom 
this point of view, a correlated initial state between PM 
and El seems to be more ‘natural’. To clarify the im¬ 
plications of the two aforementioned choices we employ 
these initial states when evaluating the BCE in setting 
(H) in the following. 


1. Factorizing initial state between PM and El 


The first initial state we consider is the one typically 
associated with factorizing initial conditions between PM 
and bath El, 


p|jj(0) = _e-/3(™^f>+i:^‘^Aata^)^ (23) 

'Zi 


where /3 = I/T is the inverse temperature and the parti¬ 
tion function Z is defined such that TrEii{/5|ij(0)} = 1. 


The second initial state we consider we call diagonal 
initial state; we thereby denote the canonical state of 
the environment EH. This state we obtain for increased 
coupling between PM and bath El, since with increasing 
coupling the equilibrium state of PM and bath will be 
given by the thermal state of the joint PM-bath environ¬ 
ment, which no longer factorizes into a PM and a bath 
part. Introducing the creation (annihilation) operators 
of the eigenmodes of the joint PM-bath system (c^), 
the global thermal state reads 

= (26) 

Here, the superscript 2) denotes a diagonal initial PM- 
bath state, implying that at t = 0 PM and bath have 
jointly evolved towards a thermal state whose occupation 
depends on the eigenenergies of the composite system. 
As before, Z is defined such that TrEii{/5gjj(0)} = I. 

For the diagonal initial state the BCF reads (r = t—t') 


a{r) = \9? E + 1)] ■ 

(27) 

Note that Eq. (27 ) is of the same form as in the stan¬ 
dard case [cf. Eq. (A5|] and can hence be written in the 
standard form Eq. ([^ with transformed SD. 
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FIG. 2. (Color online) Bath correlation functions oi{t,t')/\g\^ for different coupling strengths g and different reference times 
t'. The left column [panels (a),(c)] shows the BCF for weak PM-bath coupling, rj = 0.25, whereas the right column [panels 
(b),(d)] shows a strong coupling regime, rj = 1.0. In the first row [panels (a)-(b)] the reference time of the BCF is set to 
t' = 0, in the second row [panels (c)-(d)] t' = 32.5 A“^. Solid blue lines indicate the real part of the BCF with factorizing 
initial state [Eq. (23l], dashed green lines the real part of the BCF with diagonal initial state [Eq. (26l]. Dotted lines show the 
corresponding imaginary parts. The insets in panel (c)-(d) provide a detail of the short-time dynamics, while the inset in (a) 
shows the SD JBi(ci') (dashed red line) with the position of the PM indicated by a solid vertical line. The BCFs were calculated 
using T = 46A and = 1.5 A. 


3. Discussion 


Equations (24) and (27) allow for several observations. 
Firstly, the BCF comprises of the time evolution of the 
eigenmodes of the PM-bath environment weighted by the 
populations of the eigenmodes in the initial state. Sec¬ 
ondly, we explicitly see that in case of a diagonal initial 
state we obtain a stationary BCF, whereas in the case of 
factorizing initial conditions the BCF is non-stationary 
(for small times). This is to be expected, since for a PM- 
bath environment in thermal equilibrium the expectation 
value of any number operator (e.g. h%) should not depend 
on time — whic h is exactly what we observe if we set 
r = 0 in Eq. (27). (As a(t, t') cx 

the stationarity of (6^(t)6(t))Eii can be directly read off 
from the BCF.) 


The procedure outlined above can be generalized 
straightforwardly to, e.g., linearly-coupled chains of PMs 
of which the last PM is possibly coupled to a diagonal 
harmonic bath nans!, directly coupled PMs with inde¬ 
pendent baths [57] , or a combination of both m- As the 
BCF of the system is determined by the correlation func¬ 
tion of the PM operator directly coupled to the system, 
we simply need to adjust Heu in the above treatment; 
the calculation of the BCF then proceeds in the exact 
same manner as detailed above. 


Since the neglect of initial correlations can lead to no¬ 
ticeable differences in the dynamics [5^ UHl US] > depend¬ 


ing on the parameters of the underlying Hamiltonian, we 
now turn to the discussion of numerical examples. 


III. NUMERICAL EXAMPLES 
A. Evaluation of the transformed BCFs 


In our numerical examples, we take as spectral density 
for FI an ohmic SD with exponential cutoff (w > 0), 


>/ei(w) = rjuje 


_ —ojjA 


(28) 


where A is the cutoff frequency and ry a scaling for the 
overall coupling strength. For numerical purposes, we 
sample Jei(w) at discrete frequencies uj\. The couplings 
K\ of Eq. (14) we obtain by evaluating the quadrature 


Kx = yj (29) 

with Awa = (wa+i— a;A-i)/2 for A = 2,..., A—1; Awi = 
W 2 — wi and Awat = ujn — ujn-i- 

The sampling range is chosen such that the full SD is 
covered. For the particular cases shown, we use N = 4000 
bath oscillators for the numerical discretization, with 
equidistantly spaced frequencies, starting from 0.002 A. 
For the sake of clarity of presentation, we choose the PM 
frequency close to the maximum of the ohmic SD, setting 
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FIG. 3. (Color online) Spectral density JEii(ti;) (black, solid line) for (a) rj — 0.25 and (b) rj = 1.0 obtained from Fonrier 
transformation of the BCF (for details, see text). As the BCFs for diagonal and factorizing initial state (evaluated at tcm = 
130 A“^) coincide, they yield identical SDs. In both panels, the red, dashed line indicates the original SD of the ohmic bath, 
which has been scaled by the factor denoted in red. All other parameters are as in Fig. 


n = 1.5 A. This choice renders the coupling between PM 
and bath strongly dependent on the overall scaling of the 
SD, which is quantified by ry. 

Note that the Hamiltonian "Heii is positive for all 
parameters employed in the numerical calculations dis¬ 
cussed in this section, and that the finite recurrence time 
of the BCF is large enough to observe complete decay of 
the BCF. 


We now evaluate the BCFs Eqs. (24) and (27) with 
the described numerical procedure for different couplings 
T] and times t'. Using the SD of Eq. (28), this yields 
the BCFs displayed in Fig. There, the left column 
[panels (a),(c)] corresponds to relatively weak coupling 
77 = 0.25 whereas in the right column [panels (b),(d)] the 
PM is relatively strongly coupled to the bath El, 77 = 1.0. 
Furthermore, the first row [panels (a),(b)] show the BCF 
evaluated at = 0 while the second row [panels (c),(d)] 
show the BCF evaluated at t' = 32.5 A“^. 


As can be seen from Fig. (a) and (b), at t' = 0 
pronounced differences emerge between the two different 
initial conditions (blue vs. green) as the coupling 77 is 
increased. On the one hand, the damping of the BCF 
is increased in the presence of strong coupling, 77 = 1.0 
(note that the overall time of equilibration of the BCF 
increases as well), which results in different dynamics for 
the two different initial states. On the other hand, the 
initial values of the BCFs, a(0, 0), change, highlighting 
the increasing differences between the initial states that 
manifest themselves in the dynamics of the BCF. 

Considering the BCFs evaluated at t' = 32.5 A“^, 
shown in Fig. [2 (c) and (d), we observe that the BCFs ob¬ 
tained from different initial states look very similar, with 
the strongest difference being a transient equilibration 
dynamics present for factorizing initial conditions, which 
gets more noticeable in the strong coupling case. The 
differences found between the different initial states for 
t = t' = 0 have almost vanished for t = t' = 32.5 A“^ [cf. 


panels (a, c) and (b, d), respectively], due to the fact that 
equilibration has already taken place before t = 32.5 A“^. 

The non-stationarity of the BCF for factorizing initial 
conditions is related to what is called “initial slippage” 
if a system is coupled to a Markovian environment. In 
such systems, a non-Markovian feature can be present 
at small times due to the fact that Markovian dynamics 
for the total system requires correlations between sys¬ 
tem and environment that are not present initially if fac¬ 
torizing initial conditions are employed [UlSn]. Hence, 
slippage of initial conditions can remedy non-Markovian 
dynamics introduced by an initially uncorrelated system- 
environment state. In the same manner, the BCFs dis¬ 
played in Fig. need some time to equilibrate for fac¬ 
torizing initial conditions before reaching ‘stationarity’ 

m- 


B. Corresponding SDs 


The SD, as defined in Eq. (|^, is fully determined by 
the total Hamiltonian. Consequently, it should not de¬ 
pend on the initial state of the environment. Eor a diago¬ 
nal environment, however, the SD can be extracted from 
the BCE owing to the relation ([^ between BCE and SD. 

That is, for the diagonal initial state, the SD can 


be obtained by Eourier transforming Eq. (27) with re¬ 
spect to the time difference r and dividing by the factor 
27r(n(a;) -I- 1) [cf. Eq. (A9)]. For factorizing initial state, 
we can rewrite Eq. (24) as a function of the center of 
mass coordinate tcm = (t + ^ 0/2 and the time difference 
T = t — t' and perform a Eourier transformation with 
respect to r. 

The spectral densities corresponding to the parame¬ 
ters of Fig. are shown in Fig. For large times tcm, 
4m ^ 120 A~, the difference between the Fourier trans¬ 
formations obtained from the BCFs using a diagonal and 
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a factorizing initial state respectively, vanish. For this 
reason, we show only a single SD in Fig. 

The SDs obtained from the BCFs perfectly agree with 
the analytical result of Ref. [SI which has been obtained 
by direct transformation of the SD. (Note that in Ref. 0 
another convention for the SD has been used, i.e., the 
SD is defined as uP‘J{u}) in our convention. For a discus¬ 
sion of the advantage of the convention employed in this 
paper, see Ref. I?7ll 

As shown in Fig. for weak coupling j] = 0.25 the 
SD exhibits only a single peak at approximately the PM 
frequency fl, with a width proportional to the coupling 
77 . For large coupling 77 = 1.0, the single peak is split into 
two and the coupling strength at the PM frequency is re¬ 
duced. This illustrates that for large PM-bath coupling, 
the bath properties indeed become essential for a system 
coupled to the PM. 


C. Corresponding system dynamics 


We now turn to analyzing the effect of the features seen 
in Fig. on a system observable. To that end we specify 
the system operator L (L^) in Eq. (11) as the annihilation 
(creation) operator d (rfl) of a harmonic oscillator with 
frequency flgys, with the associated system Hamiltonian 
reading 


■Hsii = ^^sysd^d. 


(30) 


This particular set of non-Hermitian coupling operators 
allows us to conveniently evaluate the dynamics of the 
total system via diagonalization, as outlined in Sec. [IIP 


The corresponding correlation functions and 

a 2 {t,t') defined in Eqs. (j^ and ( 8 b) can be easily read 
off from Eqs. (24) and ( |27[ ). 

Setting in the above parameters 77 = 1.0, Hgys = 0.46 A, 
and requiring 7rsys(0) = (d^(0)d(0))Eii = 0, we obtain 
Fig. for two different system-PM couplings, g/A = 0.3 
and g/A = 0.08. Note that only is shown in 

Fig. since for the parameters chosen, is in¬ 

distinguishable from ai{t,t') on the scale of the figure. 
The reason for choosing flgys relatively small is that the 
steady state value of nsys{t) decreases with increasing 
Hsys, such that absolute differences in 7isys(t) are sup¬ 
pressed for large system frequencies. Likewise, we have 
to choose 77 large in order to assure that the system dy¬ 
namics is affected by the bath, being mediated via the 
PM. 


For strong coupling of the PM to both system [g = 
0.3 A, Fig. 0 (a)] and bath, we observe a marked differ¬ 
ence in the mean occupation number nsys{t) between the 
results attained from using different initial states. This 
difference highlights that in case of a strongly coupled 
PM the initial state of the environment becomes impor¬ 
tant for the transient system dynamics. In contrast, the 
equilibrium values of both nsysit) and npM(t) are inde¬ 
pendent of the initial state [3H]. If we decrease the cou¬ 
pling g to the system mode (or, similarly, the coupling 
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FIG. 4. (Color online) Bath correlation function ai(t,Q)/\g\^ 
(solid lines) and system mean occupation number nsysit) = 
id\t)d{t))-Eii (dashed lines) for different initial conditions 
and different coupling strengths g. In the upper panel (a), 
g = 0.3 A, whereas in the lower panel (b), g = 0.08 A. 
Thin blue lines indicate the BCF with factorizing initial state 
[Eq. ( |23| )], thick green lines indicate the BCF with diago¬ 
nal initial state [Eq. ( |26| )]. Only the real parts of the BCFs 
are shown. The inset shows the PM mean occupation num¬ 
ber npmit) = (b^(t)b(t)) Eli for g = 0.3 A (dashed line) and 
g = 0.08 A (solid line) for both initial states. Except for 
77 = 1.0 and Hsys = 0.46 A, the parameters of Fig.j^have been 
used. Note that ai(t, 0) is approximately Q(t, 0)/2 in Fig. 


77 to the bath oscillators), the difference in the dynam¬ 
ics decreases [cf. Fig. (b)] and the appropriate choice 
of initial conditions of the environment becomes less im¬ 
portant. The same applies for choosing flsys large, as 
this results in a lower steady-state value of Usysit) and 
consequently smaller overall deviations. 

Besides, the time it takes for the system to equilibrate 
increases as g is decreased [cf. Fig. |^(b)], as the equili¬ 
bration of the system only proceeds via interaction with 
the PM. Conversely, the PM equilibrates faster for small 
system-PM coupling g (cf. inset), since in this case the 
strong PM-bath coupling dominates the equilibration dy¬ 
namics of the PM. This again illustrates that for small 
system-PM coupling, it is indeed valid to assume a fac¬ 
torized initial environment state since the differences in¬ 
duced by the BCF rapidly vanish from the system’s point 
of view. 

For low bath temperature, the differences in the dy¬ 
namics persist, however, they become hardly observable 
due to the steady-state values (as well as the initial con¬ 
ditions) being significantly smaller as compared to high 
bath temperature. Hence, at low temperature, absolute 
deviations are reduced while relative deviations are pre¬ 
served. 

Our numerical simulations show for different initial en¬ 
vironment states pronounced differences in the transient 
dynamics of a system that via a PM strongly couples to 
an ohmic bath. Thus, any scheme sensitive to the tran- 
























sient behavior of the BCF crucially depends on the choice Appendix A: Microscopical definitions of BCF and 
of initial conditions of the total system. SD 


IV. CONCLUSIONS 


We have analytically and numerically studied the the 
BCF resulting from effectively treating a harmonic oscil¬ 
lator (PM) linearly coupled to a harmonic bath as part 
of the bath, for a factorizing and a correlated initial state 
between PM and bath respectively. We outlined the pro¬ 
cedure to analytically derive the transformed BCF and 
discussed concrete examples for regimes in which the dif¬ 
ferences in the BCFs arising from different initial states 
manifest themselves in a different dynamics of the sys¬ 
tem, which we take to be a harmonic oscillator linearly 
coupled to the PM. This establishes a simple framework 
to evaluate transformed BCFs of PMs coupled to a har¬ 
monic bath. 

In particular, we find that in the case of a correlated 
(diagonal) initial state the BCF features all the proper¬ 
ties typically assumed for a BCF (i.e. , sta tionarity, de¬ 
tailed balance and the relation Eq. (A9| between the 
Fourier transform of the BCF and the SD). By contrast, 
for a factorizing initial state these properties do not ap¬ 
ply, owed to the non-stationarity of the BCF in this case. 
Only after a transient equilibration dynamics that can 
induce different system dynamics, they are recovered. 

The significance of the BCF lies in the fact that it 
quantifies the effect of the environment — including en¬ 
vironment correlations and temperature — on the system 
degrees of freedom. Hence, our analysis complements the 
investigation of bath transformations (e.g. the mapping 
of a structured bath consisting of many harmonic oscilla¬ 
tors to a linear chain of oscillators [IIHIS]) which focused 
on the SD. 

Our findings highlight that (i) the differences between 
the BCFs for the different initial states chosen (factor¬ 
izing and diagonal initial state) are negligible for small 
PM-bath coupling or PM-system coupling respectively, 
as expected, yet (ii) these differences can have strong 
impact on the system dynamics if the PM is strongly 
coupled to both system and bath. Therefore, in case of 
a strongly coupled PM, an appropriate initial state of 
the environment has to be used for the transformation 
of the BCF when considering a system’s dynamics in the 
presence of finite temperature. This emphasizes the rel¬ 
evance of accounting for correlations in strongly coupled 
systems, which is not specific to our particular system 
[HIMIISI]- The question which initial state is to be con¬ 
sidered as appropriate cannot be answered a priori, but 
has to be answered in consideration of the specific case. 


ACKNOWLEDGMENTS 


We consider the same total Hamiltonian as introduced 
in Eq. Q, 


■Htot = +'Hs-E + ^E, (Al) 

consisting of a system part 'Hs, an environment part 'He 
[E q. ([^], and a coupling part Hs-e [Eq. ([^] accounting 
for the interaction between system and environment. 

Following Sec. EB the BCF for a Hermitian system 
coupling operator Ls is given by [cf. Eq. ([^] 

= (X! + +H.c.)\ 

/ E 

(A2) 

where the second equality sign holds if the initial environ¬ 
ment density operator /3e(0) is taken to be a function of 
the number operators cj)c^. Here, the angle brackets de¬ 
note the trace over the environment, (•)e = TrE{-pE(0)}. 


Thermal environment 

In the following, we take the environment to be in a 
thermal state. 


PE(0) = ^e-^^.^-=^-, (A3) 

where Z is defined such that TrE{/3E(0)} = 1 and /3 is 
the inverse temperature /? = I/T. 

The time evolution of the environment operators 
can be calculated by means of the Heisenberg equations 
of motion (cf. Appendix [C|), 

dtCf,(t) = (A4) 


where the time dependence refers to the interaction pic¬ 
ture with respect to He, c^(i) = 

Using the usual commutation relation [c^,cj),] = <5^^' 
and evaluating the trace in Eq. (A2) in the number basis 
I ), we find 


+ e*--(*-‘')n(w^)), (A5) 


with the mean occupation number n(a;) of the environ¬ 
ment oscillator with frequency ui defined as 


n(iij) = 


1 

- 1) E. 


We thank Gerhard Ritschel for useful discussions. 


(A6) 
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Performing a Fourier transformation with respect to 
T = t — t', we get 

Q!(a;) = 27r^ [{n{Co^) + 1) 6{uj - uj^) 




(A7) 


By means of the relation n{—ui) = — [n(a;) + 1] and the 
definition 




(A 8 ) 


we can rewrite Eq. (A7), reading 

a{uj) = 27r {n{uj) + 1) [j{uj) - ji-uj)]. (A9) 

Following Refs. [57| and [311 now define the spectral 
density J(w) as 


J{uj) = j{oj) - j{-oj) (AlO) 


and obtain, after performing the inverse Fourier trans¬ 
form and rearranging using 1 -|- 2n(a;) = coth[a;/(2T)], 
the standard expression [cf. Eq. @] 


, , du; , , 

a('r) = J ^ e a(w) 

/ OO 

duj J{uj) [n(a;) -I- 1] 

-OO 

= J dujJ{uj) ^coth cos(a;r) — *sin(wT)^ . 

(All) 


This is the well-known result for a linearly coupled har¬ 
monic environment in thermal equilibrium. 


Appendix B: Environment transformation 


Here we show how the environment Hamiltonian "Heii 
comprising both PM and bath is diagonalized. That is, 
following Ref. O we rewrite Eq. ([T^ as 


T^eii = a) Ma^ 


(HI) 


where the vector a contains all environment annihilation 
operators, 

a = ( 6 , 01 , 02 ,... (B2) 


and the matrix M all environment couplings and ener¬ 
gies. 


/H kI k.2 ••• 

Ki wi 0 • • • 0 


M = 


K2 0 UJ2 


: : 0 
yAtiv 0 • • • 0 Wat j 


(B3) 


The Hermitian matrix M can be diagonalized by means 
of a unitary transformation, 

M = SDS\ (B4) 

where the diagonal matrix D contains the eigenenergies 
of the composite bath, 



With these definitions, the new annihilation opera¬ 
tors of the environment become c = S'^d where c = 
(cq, Cl,..., cat)^. The initial creation and annihilation 
operators are obtained from the new ones via the inverse 
transformation d = Sc. Note that for a discrete number 
N of a\ operators there are N + 1 C/j, operators. 


Appendix C: Alternative derivation of the BCF on 
basis of Heisenberg equations of motion 


In this section we review a method alternative to the 
one introduced in Sec. Ill PI to calcul ate th e BCF for the 
two initial states introduced in Sec. HE [Eqs. (231 and 
(26)] by calculating the time dependence of the PM op¬ 
erator b not via diagonalization, but rather by means of 
the Heisenberg equations of motion. 

The Heisenberg equations of motion for the PM opera¬ 
tor 6 can be easily derived by evaluating the time deriva¬ 
tive 

dtb{t) = 

=-inb{t) - (Cl) 


where the time dependence as before refers to the interac¬ 
tion picture with respect to ’Heii, 6(t) = 
and 6(0) = 6 . Likewise, for a\ we have 

dtaxit) = 

= -iu}\ax{t) - iKxb{t). (C2) 


Formally integrating Eq. (C2), we find 

rt 


ax{t) = e-“^‘aA(0) - ikx f ds e-“^(‘-") 6 (s), (C3) 

Jo 

which we can insert into Eq. yielding 

dtb(t) =—iQ.b{t) — f ds W |AtApe“*‘^^^*“®^ 6 (s) 

Jo 




.* —iuJxt 


ax{0). 


(C4) 


Defining Ar(t — s) = IKApe we can write the 

solution of Eq. ( |C4| ) as 

rt 


6 (t) = ?7(t)6(0) — i y' ka®a(0) f dsU{t — s)e 
, Jo 


(C5) 
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where U(t) is determined from the integro-differential 
equation [33] 


dtU{t) =—iVtU{t) — / ds K{t — s)U{s). 


(C6) 


Note that U{t) is only defined for t > 0, with the ini¬ 
tial condition reading C/(0) = 1. Albeit for a continuous 
bath spectrum an analytic calculation of U(t) for specific 
spectral densities is possible via Laplace transforms [34] , 


we focus on a numerical scheme for solving Eq. (C61 in 


the following, relying on numerical diagonalization. 

To this end, we first define auxiliary coefficients U\{t) 
whose dependence on t is via a simple exponential and 
the integration boundary only. 


Ux{t) = Kxf 
Jo 


(C7) 


By means of this definition, we are able to cast Eq. (C6) 
into a set of coupled equations, 


dtUit) = -iVlU(t) - ^ K*xUx{t), 

(C8a) 

dtUx(t) = KxUit) - iuxUxit). 

(C8b) 

Introducing the vector u{t) = {U{t),Ui{t),. 
we can rewrite Eqs. (C8a) and (C8b| as 


dtu{t) = —iGu{t), 

(C9) 


with the matrix G given by 


f Q —ini —1^2 

IKi OJi 0 


G = 


iK2 0 


UJ2 


-iK%\ 

0 


\iK]s[ 0 


0 

0 Wat / 


(CIO) 


The Hermitian matrix G can be diagonalized via the 
transformation GT = D, yielding the eigenvalues 
Defining u = Tv, the differential equation Eq. (C9) be¬ 
comes 


(Cll) 


which has the simple solution v{t) = with 

u(0) = Tfu(O), where u(0) = (1,0,..., 0)^. Tak¬ 
ing the first component of the full expression u{t) = 
Te“*^*r'^u(0), which gives U{t), and writing [■u(O)]^ = 
So/_i, we arrive at the simple expression 

C/(t) = (C12) 


The above reformulation allows one to map the so¬ 
lution of an integro-differential equation onto an eigen¬ 
value problem which involves only matrix multiplication 
and diagonalization, which is numerically more robust 
with respect to the numerical time step than straightfor¬ 
ward numerical integration of Eq. ( |C6| ). Note that the 
approach used requires K{t — s) to be given as a sum 
of exponentials that reproduce when being differentiated 
with respect to time. 


We can now employ the above results to obtain alter¬ 
native analytic expressions for the BCFs of the two initial 
states introduced in Eqs. (23) and (26). 


Eor factorizing initial conditions, Eq. (23), the trace in 
Eq. (22) is readily evaluated, yielding 


a{t,t')/\g\^ = U{t)U*{t') (n(D) -f 1) -f U*{t)U{t')n{n) 


E 


K\\ 



pt r 

/ ds 

/ ds' 

Jo 

lo L 


U{t - s) 


X U*(t' - {n{ujx) + 1) 

+ U*{t- s)U{t' - . 

(C13) 


Here, n(a;) is the mean occupation number of the 
harmonic oscillator of frequency ui, as introduced in 
Eq. (|A6|. 


Eor diagonal initial conditions, Eq. (26), it is advis¬ 
able to first linearly transform the operators b, ax into 
the eigenbasis of the joint environment, outlined in Ap¬ 
pendix]^ We find 


dtv{t) = —iDv{t), 
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a(t,i')/l5l" = E [U{t)U*{t') {n{Q^) + 1) + C/*(t)C/(t')ri(cD^)] + E f ds f ds' 

M A.r,M "'O Jo 




(n{ujfj_) + 1) + U* {t — s)U{t' — s')e 






X.fi 

A,// 


C/(i) [ dsU*{t' - s)e^^^^n{ujf,) + 1) + U{t') f dsU* it - s)e'‘^^‘‘n{uj^) 

Jo Jo 

U*{t') [ dsUit-s)e-^‘^^^{n{CJf,) + l) + U*{t) f dsU(t' - s)e-'-‘^^^n{Cb^) 


■ (C14) 


Note that the above results assume a Hermitian system 
operator L of the system-PM coupling [cf. Eq. 0]- 

The procedure allowing us to arrive at Eqs. ( |C13| ) and 
(C14) is straightforward: We (i) derived the Heisenberg 


equations of motion for the PM operator, (ii) linearly 
transformed the time-independent {t = 0) operators of 
the Heisenberg equations of motion into the operators 
with respect to which the initial state is diagonal (cf. 
Appendix H , and (iii) evaluated the trace of the BCE 
[Eq. (22)] in this basis. 


Analogously to the procedure described in Sec. HD 


the scheme summarized above can be generalized to, e.g., 
linearly-coupled chains of PMs with the last PM(s) pos¬ 
sibly coupled to a terminating bath nzms], direct cou¬ 
pling of PMs to independent baths or a combina¬ 
tion of both m- Since the terminating baths are usu¬ 
ally taken to be independent, the BCFs of the PMs can 
be derived by successively solving the Heisenberg equa¬ 
tions of motion following the above approach. Starting at 
the PM which is coupled to the system and subsequently 
transforming the PM and bath operators into the basis 
in which the initial states are diagonal, the BCFs can be 
evaluated with respect to the given initial states. 
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